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ABSTRACT 


The  NASTRAN  structural  analysis  program  has 
been  modified  so  that  its  TRAPAX  and  TR1AAX  finite 
elements  may  handle  piezoelectric  material 
properties.  These  elements  are  solid,  axisymmetric 
ring  elements  capable  of  handling  non-axisymmetric 
loads.  The  natural  frequencies  of  a  piezoelectric 
hollow  cylinder  and  a  piezoelectric  disk  with 
electroded  surfaces  were  computed  with  this 
modified  version  of  NASTRAN  and  were  found  to  be 
in  excellent  agreement  with  experimental  and 
other  numerical  results. 


ADMINISTRATIVE  INFORMATION 

This  project  was  funded  by  the  Naval  Underwater  Systems  Center,  New 
London  Laboratory,  under  the  direction  of  Mr.  A.D.  Carlson.  The  work  was 
performed  under  Task  Area  RR0310301  and  Work  Unit  1844-130. 


INTRODUCTION 


The  analysis  of  sonar  transducers  requires  accounting  for  the 

effects  of  their  piezoelectric  materials.  The  theory  for  these  materials 

1  2* 

couples  the  structural  displacements  to  the  electric  potentials.  ’ 

This  theory  has  heen  incorporated  into  the  TRAPAX  and  TRIAAX  finite 

3 

elements  of  the  NASTRAN  structural  analysis  computer  program.  These 
elements,  trapezoidal  and  triangular  in  cross-section  respectively,  are 
solid,  axisymmetric  rings  whose  degrees-of-f reedom  are  expanded  into 
Fourier  series,  thus  allowing  nonaxisymmetr ic  loads. 

The  complete  details  of  the  TRAPAX  and  TRIAAX  finite  elements  and  of 

12  3 

finite  element  piezoelectric  theory  are  well-documented  ’  ’  and  will  not 
be  repeated  here.  However,  some  of  the  ideas  will  be  briefly  discussed. 
The  constitutive  relations  of  the  material  may  be  written  as 

E, 


{o}  He 
(D}|  |_[e  ] 


[e] 

-U-S] 


(1) 


*A  complete  listing  of  references  is  given  on  page  25. 
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where  {a}  =  stress  components  *  la  ,  a  ,  a  ,  a  „ ,  a  ,  F 

L rr  zz  06  rz’  r6  z6J 

l  T 

{Dl  =  components  of  electric  flux  density  =  Id  ,  D  ,  D  ,  \ 

J  *-  rr  zz  0 0-i 

{e}  =  mechanical  strain  components 

{ E }  =  electric  field  components 

£ 

fc  ]  =  elastic  stiffness  tensor  evaluated  at  constant  electric  field 
[ e ]  =  piezoelectric  tensor 

S 

[e  ]  =  dielectric  tensor  evaluated  at  constant  mechanical  strain 


The  displacement  vector  of  a  point  within  an  element  is  taken  to  be 


where  u,  v,  and  w  are  the  ring  displacements  in  the  radial,  tangential, 

and  axial  directions,  respectively,  and  <p  is  the  electric  potential.  The 

latter  degree-of- f reedom  is  taken  to  be  the  fourth  degree-of-f reedom  at 

each  ring.  Each  of  these  quantities  is  expanded  into  a  Fourier  series 

with  respect  to  the  azimuth  position  0.  The  Fourier  series  for  the 

electric  potential  f  has  the  same  form  as  the  Fourier  series  for  radial 

3 

displacement  u,  as  given  in  the  NASTRAN  Theoretical  Manual. 

The  "stiffness"  matrix  for  the  harmonic  is 


[K(N)]  -  7T  /  j  [B(N)  JT  [c  ]  [e]  [B<N)]  rdrdz  (3) 

r  z  [e}L  -[e  ] 

U.  J 

(N) 

where  [B  ]  is  the  matrix  of  "strain"-"displacement"  coefficients  for 
the  harmonic. 

Equations  (2)  and  (3)  indicate  that  the  matrix  equation  to  be  solved 
for  static  analysis  may  be  partitioned  as  follows: 
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where 

{5}  - 

Lv  v  v  ••••  V 

lT 

v  ,  w  J 
n  nJ 

U)  = 

T 

LV  •••» 

{V  - 

vector  of  structural 

forces 

and 

{V  - 

vector  of  electrical 

charges 

Note,  however,  that  the  program  assumes  that  the  electric  potential 
is  the  fourth  degree-of-f reedom  of  grid  point  i. 

Both  lumped  and  consistent  mass  matrices  are  available  and  are  of 
standard  structural  form,  i.e.,  the  mass  matrix  does  not  couple  the 
structural  and  electrical  unknowns.  The  structural  damping  matrix  also 
is  of  standard  structural  form.  Both  point  charges  and  surface  charges 
are  also  available. 


NASTRAN  MODIFICATIONS 

Numerous  modifications  were  required  to  incorporate  the  piezoelectric 
material  capability  into  NASTRAN.  These  modifications  included  new  bulk 
data  cards,  new  outlooks  on  existing  data  cards,  and  subroutine  changes, 
some  of  which  were  the  result  of  existing  errors.  All  modifications 
were  made  to  NASTRAN  Level  17.5. 

DATA  CARDS 

Four  new  bulk  data  cards  have  been  introduced  into  NASTRAN,  and  two 
existing  bulk  data  cards  can  be  viewed  in  a  different  light  to  accommodate 
piezoelectric  materials.  Also,  the  NASTRAN  card  contains  a  new  parameter. 

New  Cards 

The  four  new  bulk  data  cards  adder  to  the  program  describe 
piezoelectric  material  properties.  These  properties  are  frequently 
described  by  the  following  matrices: 
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[SE]  = 


where  S 
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(5) 


(6) 


(7) 


The  matrices  in  Equation  (1)  are  computed  as  follows: 

[cE]  =  [SE]-1 
[e]  =  [d][cE] 


(8) 

(9) 


and  [c  ]  is  given  by  Equation  (7) . 

Two  of  the  new  data  cards,  MATPZ1  and  MATPZ2,  describe  the 

piezoelectric  material  properties  in  two  different  ways.  MATPZ1  is  used 

to  specify  the  parameters  in  Equations  (5)-(7).  MATPZ2  is  more  general 

E  S 

and  allows  the  user  to  enter  the  full  matrices  [c  ],  [e],  and  [t  1.  The 
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E  S 

only  assumption  concerning  these  matrices  is  that  [c  ]  and  [e  ]  are 
symmetric.  Figures  1  and  2  describe  MATPZ1  and  MATPZ2. 

CAUTION:  Piezoelectric  electric  material  properties  are  usually 

specified  with  respect  to  a  standard  set  of  material  axes  1,  2,  3.  For 
axisymmetric  solids,  direction  1  coincides  with  the  z-axis  and  direction  2 
coincides  with  the  0-axis.  Polarization  direction  3  may  vary  in  the 
R-Z  plane  and,  for  radial  polarization,  coincides  with  the  R-axis .  When  a 
user  specifies  properties  on  a  MATPZ1  card,  the  transformation  from  the 
1,  2,  3  directions  to  the  R,  Z,  0  directions  is  performed  by  NASTRAN. 
However,  such  a  transformation  is  not  performed  by  NASTRAN  when  the 
MATPZ2  card  is  used.  Also,  the  ordering  of  the  components  of  the  stress 
and  strain  vectors  is  somewhat  different  for  conventional  piezoelectric 
specifications  and  for  NASTRAN.  The  difference  is  that  the  ordering  of 
the  Z-0  and  R-Z  shears  is  interchanged.  Once  again,  NASTRAN  performs  the 
transformation  for  MATPZ1,  but  not  for  MATPZ2. 

The  other  new  data  cards,  MTTPZ1  and  MTTPZ2 ,  allow  the  values  on 
MATPZ1  and  MATPZ2  to  be  temperature-dependent.  (However,  the  TRAPAX 
and  TR1AAX  elements  have  not  yet  been  modified  to  handle  the  combination 
of  thermal  loads  and  piezoelectric  materials.)  Figures  3  and  4  describe 
MTTPZ1  and  MTTPZ2 . 

Existing  Cards 

Point  and  surface  charges  may  be  specified.  These  charges  are 
analogous  to  structural  point  loads  and  pressures,  respectively,  and  are 

entered  into  (F  }  in  Equation  (4).  Since  the  electric  potential  is 

4> 

associated  with  degree-of-f reedom  4,  point  charges  may  be  applied  to 
specific  harmonics  with  MOMAX  bulk  data  cards  or  may  be  specified  by 
MOMENT,  MOMENT 1 ,  or  MOMENT 2  cards  applied  to  POINTAX  points.  In  the 
latter  case,  the  direction  of  the  "moment"  must  be  about  the  radial 
direction,  i.e.,  degree-of-f reedom  4.  The  MKS  unit  of  the  point  charge 
is  coulombs. 

The  PRESAX  bulk  data  card  is  used  to  specify  surface  charges.  How¬ 
ever,  in  order  to  distinguish  between  surface  charges  and  structural 
pressure  loads  within  the  same  problem,  the  first-specified  ring 
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Input  Data  Card  MATPZ1  Material  Property  Definition 


Description:  Defines  the  material  properties  for  linear,  temperature- 

independent,  piezoelectric  materials. 

Format : 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

MATPZL 

MID 

£ 

*33 

E 

S44 

sE 

12 

,  S13  j 

d31 

+.4 

1  +A 

H 

d15 

S  / 

c33/eO 

RHO 

A 

TREF 

GE 

Field 

MID 


LI 


-  d 


S  /■- 

'll  'O’ 


15 

S  h 
33'  '0 


Contents 

Material  identification  number  (Integer  >  0) 

12 

Piezoelectric  constants  multiplied  by  10  (Real) 

Piezoelectric  constants,  where  £„  is  taken  to  be 

-12  U 
8.854  x  10  farad/meter  (Real) 


RHO 

A 

TREF 

GE 


Mass  density  (Real) 

Thermal  expansion  coefficient  (Real) 

Thermal  expansion  reference  temperature  (Real) 
Structural  element  damping  coefficient  (Real) 


Remarks : 

1.  MID  must  be  unique  with  respect  to  all  other  mates  ’"al  cards. 

2.  MATPZ1  materials  may  be  made  temperature-dependent  by  use  of 
the  MTTPZ1  card. 

3.  MATPZ1  may  be  referenced  only  by  PTRAPAX  and  PTRIAAX  cards. 

E 

4.  Matrix  [S  ]  must  be  nonsingular. 


Figure  1  -  MATPZ1  Bulk  Data  Card 
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Input  Data  Card  MATPZ2  Material  Property  Definition 


Description:  Defines  the  material  properties  for  linear,  temperature- 

independent,  piezoelectric  materials. 

Format : 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

MATPZ2 

MID 

E 

cu 

E 

C12 

E 

C13 

F. 

C14 

E 

C 1 3 

E 

C16 

E 

C22 

-A 

+A 

E 

C23 

E 

C24 

n 

E 

C26 

E 

C33 

E 

C  34 

E 

C  35 

E 

C  36 

+B 

+B 

E 

C44 

E 

C45 

•4 

-cl 

55 

E 

C3b 

.E 

1  66 

°11 

C12__ 

+C 

+C 

e13 

e14 

e15 

e16 

e21 

e22 

e23 

e24 

+D 

+D 

e25 

e26 

e31 

e32 

0 

33 

e34 

e35 

e36 

+E 

+E 

S 

£11 

S 

C12 

S 

6 13 

S 

£22 

s 

£23 

S 

E  33 

RHO 

AX 

+F 

+F 

AY 

TREF 

CE 

Field 

MTD 

E  S 
CU  ~  '  3  3 

RHO 

AX, AY,AZ 
TREF 
CE 


Contents 

Material  identification  number  (Integer  '•  0) 
Piezoelectric  constants  (Real) 

Mass  density  (Real) 

Thermal  expansion  coefficients  (Real) 

Thermal  expansion  reference  temperature  (Real) 
Structural  element  damping  coefficient  (Real) 


Remarks : 

1.  MID  must  be  unique  with  respect  to  all  other  material  cards. 

2.  MATP/.2  materials  may  be  made  temperature-dependent  by  use  ol  tin 
MTTPZ2  card. 

3.  MATPZ2  may  be  referenced  only  by  PTRAPAX  and  PTR1AAX  cards. 

4.  Sec  CAUTION  on  page  5. 


Figure  2  -  MATPZ2  Bulk  Data  Card 
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Input  Data  Card  MTTPZ1  Material  Temperature  Dependence 

Description:  Specifies  table  references  for  those  piezoelectric  material 
properties  on  a  MATPZ1  card  which  are  temperature-dependent. 

Format : 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

MTTPZ1 

MID 

Rl 

R2 

R3 

R4 

R5 

R6 

R7 

+A 

+A 

R8 

R9 

RIO 

Rll 

R12 

R13 

R14 

Field  Contents 

MID  Material  property  identification  number  which 

matches  the  identification  number  on  some  basic 
MATPZ1  card  (Integer  >  0) 

Ri  References  to  table  identification  numbers  for  the 

corresponding  fields  on  the  MATPZ1  card 
(Integer  i  0  or  blank) 


Remarks : 


1. 

quantity 

2. 

3. 

values . 
then  two 


Blank  or  zero  entries  mean  no  Cable  dependence  of  the  referenced 
on  the  basic  MATPZ1  card,  and  the  quantity  remains  constant. 
TABLEM1 ,  TABLEM2 ,  TABLEM3 ,  and  TABLEM4  type  tables  may  be  used. 
Material  properties  given  on  the  basic  MATPZ1  card  are  initial 
If  two  or  more  quantities  are  to  retain  a  fixed  relationship, 
or  more  tables  must  be  input  to  define  the  relationship. 


Figure  3  -  MTTPZ1  Bulk  Data  Card 


Input  Data  Card  MTTPZ2  Material  Temperature  Dependence 

Description;  Specifies  table  references  for  those  piezoelectric  material 
properties  on  a  MATPZ2  card  which  are  temperature-dependent. 

Format: 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

MTTPZ2 

MID 

Rl 

R2 

R3 

R4 

R5 

R6 

R7 

+A 

(etc.) 

+F 

— 

R48 

R49 

R50 

— 

R51 

| 

Field  Contents 

MID  Material  property  identification  number  which 

matches  the  identification  number  on  some  basic 
MATPZ2  card  (Integer  >  0) 


Ri  References  to  table  identification  numbers  for  the 

corresponding  fields  on  the  MATPZ2  card 
(Integer  ^  0  or  blank) 


Remarks : 


1. 

quantity 

2. 

3. 

values . 
then  two 


Blank  or  zero  entries  mean  no  table  dependence  of  the  referenced 
on  the  basic  MATPZ2  card,  and  the  quantity  remains  constant. 
TABLEM1 ,  TABLEM2 ,  TABLEM3 ,  and  TABLEM4  type  tables  may  be  used. 
Material  properties  given  on  the  basic  MATPZ2  card  are  initial 
If  two  or  more  quantities  are  to  retain  a  fixed  relationship, 
or  more  tables  must  be  input  to  define  the  relationship. 


Figure  4  -  MTTPZ2  Bulk  Data  Card 
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identification  number  on  the  PRESAX  card  (field  4)  must  be  made  negative 
if  a  surface  charge  is  desired.  Modifications  have  been  made  to  NASTRAN 
to  allow  a  negative  ring  identification  number  only  when  the  parameter 
SYSTEM(78)  is  set  to  1  on  the  NASTRAN  card. 

The  NASTRAN  card  allows  the  user  to  override  various  NASTRAN  system 
parameters.  SYSTEM(78),  set  aside  for  specifying  piezoelectric  materials, 
has  a  default  value  of  zero,  which  implies  that  no  piezoelectric 
materials  are  allowed.  If  SYSTEM(78)  =  1,  piezoelectric  materials  are 
allowed  and  coupling  occurs  between  the  structural  and  electric  degrees- 
of-freedom.  If  SYSTEM(78)  =  2,  piezoelectric  materials  are  allowed,  but 
no  coupling  occurs. 

Setting  SYSTEM(78)  to  its  proper  value  is  important  for  a  number  of 
reasons : 

1.  If  SYSTEM(78)  =  0,  no  piezoelectric  materials  are  expected, 
and  MATPZ1  and  MATPZ2  cards  will  not  be  searched. 

2.  If  SYSTEM(78)  j*  1 ,  a  negr*-ive  ring  identification  number  is  not 
allowed  on  the  PRESAX  card. 

3.  If  SYSTEM(78)  ^  1,  NASTRAN  will  automatically  constrain  degree- 
of-freedom  4  (the  electric  potential)  at  each  ring  for  the  zero  harmonic 
in  the  AXISYMMETRIC  =  COSINE  case. 

4.  If  SYSTEM(78)  =  2,  some  time  will  be  saved  in  generating  the 
"stiffness"  matrix  over  the  time  for  the  SYSTEM(78)  =  1  case. 

5.  If  SYSTEM(78)  4  1,  degrees-of-f reedom  4,  5,  and  6  must  be 
removed  from  the  problem  via  the  SPCAX  or  RINGAX  cards.  If  SYSTEM(78)  =  1 
only  degrees  of  freedom  5  and  6  must  be  removed. 

SUBROUTINES 

Numerous  modifications  have  been  made  to  NASTRAN  Level  17.5  to  allow 
piezoelectric  materials,  point  and  surface  charges,  and  complex  stress 
and  complex  force  output  for  TRAPAX  and  TRIAAX  elements,  and  to  correct 
existing  errors  related  to  these  elements.  Tne  modifications  for  each 
subroutine  are  described  in  this  section. 
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DELKLS 

The  modification  to  this  subroutine  allows  for  finer  finite  element 
meshes  with  the  TRIAAX  and  TRAPAX  elements  by  tightening  a  check  on  the 
difference  between  radial  coordinates  of  the  grid  points  of  an  element. 

EMA 

If  a  structural  damping  coefficient  is  present,  only  the  uncoupled 
structural  terms  of  the  stiffness  matrix  are  multiplied  by  it.  All  the 
piezoelectric  terms  of  the  stiffness  matrix  are  set  to  0. 

GPTABD 

An  array  C0MP3  was  added  to  describe  the  manner  in  which  the  real 
and  imaginary  parts  of  complex  stresses  and  forces  are  related  to  eacli 
other  for  output  for  the  TRIAAX  and  TRAPAX  elements.  Also,  the  E  array 
for  these  elements  was  updated  to  point  to  the  new  information  in  COMPJ. 
GP3A 

Provision  was  made  to  allow  a  negative  ring  identification  number 
on  the  PRESAX  card  to  indicate  a  surface  charge  rather  than  a  structural 
pressure. 

IFP 

Since  the  MATPZ2  and  MTTPZ.2  bulk  data  cards  are  fixed  ended  and 
contain  more  than  50  items  each,  dimensions  for  arrays  I,  M,  MF ,  and  Ml 
in  AFPDTA/were  increased  to  100.  Also,  standard  modifications  required  to 
access  new  bulk  data  cards  were  made. 

IFPDCO 

As  in  IFP,  dimensions  for  /IFPDTA/  were  increased. 

IFP  3 

Modifications  to  IFP3  were  required  to  allow  negative  ring  identifi¬ 
cation  numbers  on  PRESAX  cards  for  surface  charge  density  loading. 

IJFP3B 

In  standard  structural  axisymmetric  problems,  IFP3B  constrains 
degrees-of-freedom  2,  4,  and  6  in  the  zeroth  harmonic  for  the  COSINE  case. 
Since  voltage  is  degree-of-f reedom  4,  IFP3B  was  modified  to  constrain 
only  degrees-of-freedom  2  and  6  in  piezoelectric  problems. 
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IFS1P 

Dimensions  in  /IFPDTA/  were  increased,  and  checks  were  adaed  for 
the  data  on  new  cards  MATPZ1,  MATPZ2,  MTTPZ1,  and  MTTPZ2 . 

IFS2P 

Dimensions  in  /IPDTA/  were  increased  as  in  IFP. 

IFS3P 

Dimensions  in  /IFPDTA/  were  increased  as  in  IFP,  and  changes  were 
made  to  allow  the  PRESAX  card  to  have  a  negative  ring  identification 
number  in  piezoelectric  problems. 

IFS4P,  1FS5P,  IFXDBD 

Dimensions  in  /IFPDTA/  were  increased  as  in  IFP. 

IFXOBD 

In  IFX1BD,  all  available  spaces  for  bulk  data  cards  were  used  by 
the  standard  program.  In  order  to  increase  the  number  of  spaces 
available,  space  had  to  be  increased  in  the  IB  and  IC  arrays  in  IFXOBD, 
which  contains  bit  specifications  for  the  bulk  data,  parameter,  and  case 
control  cards  in  IFX1BD.  Associated  bit  specifications  also  had  to  be 
entered . 

IFXIBD  -  IFX7BD 

Standard  additions  were  made  for  new  bulk  data  cards  MATPZ1,  MATPZ2 , 
MTTPZ1 ,  and  MTTPZ2 . 

OFjVZZZ  (OFP) 

Two  errors  were  corrected.  Output  for  real,  sort  2,  transient 
analyses  of  axisymmetric  problems  is  now  printed  rather  titan  punched. 
Also,  if  the  problem  is  a  real,  sort  2,  statics  analysis,  the  entire 
label  of  the  heading  is  deleted  rather  than  only  part  of  it.  Modifica¬ 
tions  allow  for  the  output  from  axisymmetric  problems  of  new  displacement 
vector  headings  for  real,  sort  2,  transient  problems;  tor  real,  sort  2, 
static  problems;  and  for  complex,  sort  2,  frequency  response  problems. 
Special  packing  formats  added  to  0FP5BD  are  also  checked. 

0FP1A 

This  subroutine  contains  FORMAT  statements  for  output  headings.  Two 
existing  FORMAT  statements  were  modified  to  allow  the  output  of  charges 
and  flux  densities.  Thirteen  FORMAT  statements  were  added  to  provide  a 


complete  set  of  output  headings.  These  FORMAT  statements  include 
headings  for  stresses,  forces,  loads,  and  displacements  for  complex 
output,  sort  1  and  sort  2,  and  real,  sort  2,  transient  and  static  output. 
OITIBD 

This  BLOCK  DATA  subprogram  is  pointed  to  from  one  of  10  tables, 
depending  on  the  type  of  stresses  or  forces  being  processed  for  output, 
ft  points  into  0FP5BD  to  construct  the  correct  FORMAT  statement.  Modifi¬ 
cations  and  additions  allow  for  the  addition  of  charges  and  flux  densitie 
to  the  output  for  any  type  of  real  or  complex  problem  containing  TRAPAX 
and  TR1AAX  elements. 

0FP5BD 


Pieces  of  FORMAT  statements  from  this  block  data  were  modified  or 
added  to  allow  the  output  of  charges  and  flux  densities. 


0F2PBD,  0F1PBD,  0F3SBD,  0F4PBD,  0F5PBD, 

0F6PBD ,  0F7PBD,  0F7SBD,  0F8PBD 

These  BLOCK  DATA  subprograms  are  indirectly  pointed  to  from 
OFPZZZZ  and,  in  turn,  point  to  OFP1BD  and  0FP1  or  0FP1A.  They  modify  or 
provide  the  capability  for  output  of  the  following  types  of  problems: 
complex  stresses,  sort  1;  real  stresses,  sort  2  (not  statics);  real 
stresses,  sort  2  (statics);  complex  stresses,  sort  2;  real  forces,  sort 
complex  forces,  sort  1;  real  forces,  sort  2  (not  statics);  real  forces, 
sort  2  (statics);  and  complex  forces,  sort  2. 

I ’REMAT Z  (PREMAT) 


This  subroutine  now  allows  piezoelectric  material  properties  to  be 
read  from  the  MATPZ1,  MATPZ2 ,  MTTPZ1,  and  MTTPZ2  cards  already  described, 
it  the  piezoelectric  material  properties  are  specified  on  a  MATi’Zl  card, 
they  are  transformed  to  the  form  found  on  a  MATPZ2  card  (Equations  8 
and  0).  All  the  piezoelectric  material  constants  are  stored  in 
PZOlfT(l-51) ,  which  is  part  of  named  common  MATPZ. 

PRESAX 


A  negative  ring  identification  number  specified  on  a  PRESAX  card 
indicates  a  surface  charge  load  rather  than  a  structural  pressure  load. 
Tills  subroutine  has  been  modified  to  identify  surface  charge  loads,  and 
to  compute  the  charge  to  be  applied  at  the  rings  specified  on  the  PRESAX 
card  . 
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SDR2D 

The  pointers  into  the  displacement  vector  file  were  modified  to 
handle  any  type  of  complex  problem  for  the  TR1AAX  and  TRAPAX  elements. 
SDR2E 

Modifications  were  made  to  handle  the  complex  stress  and  force 
output  from  STRAX2  or  STPAX2.  The  imaginary  part  of  the  complex  stress 
and  force  output  is  stored  in  variables  ISAVES  and  ISAVEF.  The  dimension 
of  each  of  these  variables  was  doubled  to  handle  the  additional  charges 
and  flux  densities.  These  variables  are  also  part  of  a  new  named  common 
I SAVE. 

STPAX1 

This  subroutine  forms  the  stiffness  and  stress  matrices  for  TRAPAX 

elements  for  stress  and  force  recovery.  Two  of  the  new  matrices  to  be 

formed  are  [K.  ]  and  [K..],  the  piezoelectric  stiffne  s  matrix  partitions 
041  <p  <p 

of  Equation  (4).  [K.  .]  is  not  formed  because 

<po 

‘V  -  tKs/  <10 


To  generate  [K.r  1  and  [K,l  a  new  4x4  transformation  matrix  from  field 

o  <p  <p<p 

coordinates  to  grid  point  degrees  of  freedom  is  formed.  It  is  stored  in 

GBP  and  is  similar  to  GABABQ  used  to  generate  [K  ] .  A  piezoelectric 

E  66 

materials  matrix  (Equation  (1))  is  formed.  [c  ]  is  stored  in  EE(1-21), 
[e]  is  stored  in  EE(37-54),  and  [e^]  is  stored  in  EE(55-63)  .  [B^^]  ol 

Equation  (3)  is  never  explicitly  formed.  Rather,  [K,  )  and  [K  ]  are 

o  4‘  v 

generated  term  by  term  incorporating  the  matrix  triple  product  and 
integration  over  r  and  z  found  in  Equation  (3).  The  new  untransformed 
partitions  of  the  stiffness  matrix  are  stored  in  ACURL(145-208) .  Once 
these  partitions  have  been  transformed,  they  are  stored  in  AKUPH  and 

AKPH2,  corresponding  to  [ K _ ^ ]  and  [K  1,  respectively.  Three  new 

o<p  q> 

partitions  of  the  master  stress  matrix  must  also  be  generated.  The 
augmented  stress  matrix  is  of  the  form: 


lS0^30xl2 

[SPl]30x4 

_tSP2115xl2 

[SP3'l5x4_ 

(ID 
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where  ( S  I  =  augmented  master  stress  matrix 

4 

1  1  =  original  stress  matrix 

l^pj  I » lsp? 1 » 1SP3  1  =  piezoelectric  partitions  of  the  master  stress 

matrix 

IS  I.  lSp.,  I,  and  IS  1  are  stored  in  SELP1,  SELP2 ,  and  SELP3, 

respectively.  To  generate  the  three  new  stress  matrix  partitions,  another 

1x4  transformation  matrix  is  required.  It  is  stored  in  WJP  and  is 

similar  to  WJ  used  to  generate  [S  ].  The  new  stress  and  stiffness  matrix 

partitions  are  passed  to  STPAX2  through  a  named  common  SDR2X5.  The 

complete  NASTRAN  method  for  formulating  the  stiffness  and  stress  matrices 

4 

is  described  in  the  NASTRAN  Programmers  Manual. 

ST PAX  2 

With  the  augmented  stress  and  stiffness  matrices  from  STPAX1, 
lorces,  charges,  stresses,  and  flux  densities  are  computed  according  to 
Equation  (4)  and 


so'  lSPl 
SP2*  'SP3^ 


1  , 

1 

IS  j 

IT) 

,  - 

5 

Ift 

IT  J 

(12) 


whete  iT  ;  -  stress  (second-order  tensor) 
and  -.T  V  =  flux  density  (vector) 


■F  ,  F  'T  r,  and  {T  i  are  stored  in  EFORC,  ECHRC,  ESTRES,  and  EFLUX, 
respectively.  These  results  for  each  harmonic  and  angle  are  then 
accumulated  in  BLOCK  if  they  are  real  output  or  the  imaginary  part  of 
complex  output.  if  they  are  the  real  part  of  complex  output,  they  are 
accumulated  in  CLOCK.  However,  the  results  for  each  harmonic  at  zero 
angle  are  always  stored  in  STRES  and  FORCE,  whether  they  are  real  or 
comp  1  ex . 

Si  PAX  1 

Stresses  and  forces  from  BLOCK  and  CLOCK  are  stored  in  FORCE  and 
STRESS  or  SAVEF  and  SAVES.  In  real  problems,  output  from  BLOCK  is  trans- 
terred  to  CLOCK  and  stored  in  FORCE  and  STRESS.  In  complex  problems,  the 
imaginary  part  of  the  output  stored  in  BLOCK  is  saved  in  SAVEF  and  SAVES. 


i  5 


The  real  part  of  complex  output  from  CLOCK  is  stored  In  FORCE  and 
STRESS. 

STRAX1 

The  modifications  are  identical  to  those  in  STPAX1  except  that  the 
TRIAAX  element  has  three  grid  points  rather  than  four.  Also,  for  the 
TRAPAX  element,  stresses  are  computed  at  the  four  grid  points  and  the 
centroid,  rather  than  just  at  the  centroid  as  for  the  TRIAAX  element. 
These  considerations  change  the  size  of  the  following  matrices  to  those 
indicated:  transformation  matrix  GABABP(3x3),  ACURL(117),  [S^](6x9), 
[Sp^](6x3),  [Sp2l(3x9),  [Sp^](3x3),  and  transformation  matrix  WJP(3x3). 

STRAX2 

The  modifications  are  identical  to  those  in  STPAX2,  with  the 
exception  noted  for  STRAX1. 

STRAX3 

The  modifications  are  identical  to  those  in  STPAX3,  with  the 
exception  noted  for  STRAX1. 

TRAPAX 

Modifications  to  this  subroutine  are  identical  to  those  for  STPAX1 

for  the  generation  of  the  stiffness  matrix  partitions  [K^],  [K^],  and 

[R,,],  in  addition,  the  stiffness  matrix  returned  from  this  subroutine 

must  be  in  a  usable  form  to  compute  grid  point  displacements.  Therefore, 
T 

[Krr],  [K..J,  [ K - . ]  ,  and  [K,,]  must  be  assembled  into  a  total  element 

00  o(J>  ocp  <p<p 

stiffness  matrix.  Since  the  electric  potential  is  the  fourth  degree- 

T 

of-freedom,  [K.J,  [Kfi]  ,  and  [K  1  must  be  merged  with  [K _ - ] .  These 

0  <p  0<J)  <p(p  00 

partitions  are  inserted  in  every  fourth  column  and  row.  No  modifications 
were  required  in  the  generation  of  the  mass  matrix  to  account  for  the 
piezoelectric  capability.  However,  two  existing  errors  in  TRAPAX  were 
corrected.  In  the  first  of  these,  the  terms  of  the  lumped  mass  matrix 
were  not  divided  by  four.  Without  dividing  by  four,  the  total  lumped 
mass  was  assigned  to  each  grid  point  rather  than  only  one-fourth  of  it. 

In  the  second,  an  error  in  the  coding,  the  consistent  mass  matrix  was 
not  transformed  from  field  coordinates  to  grid  point  degrees-of-freedom. 


TRIAAX 

Modifications  to  this  subroutine  are  identical  to  those  for  STRAX1 

for  the  generation  of  the  stiffness  matrix  partitions  [K,.],  [K,  ],  and 

oo  09 

[K,t].  The  formulation  of  the  element  stiffness  matrix  and  the  mass 

r  *' 

matrices  follows  the  same  procedure  as  does  the  TRAPAX  subroutine.  One 
existing  error  was  corrected  in  TRIAAX.  The  terms  of  the  lumped  mass 
matrix  were  not  divided  by  three.  Without  dividing  by  three,  the  total 
lumped  mass  was  assigned  to  each  grid  point  rather  than  only  one-third 
of  it . 

TTl.PGE 

Modifications  update  the  NASTRAN  title  page  to  Level  17.5.1  from 
Level  17.5.0  and  to  a  system  generation  date  of  8/15/79  from  11/30/78. 

OVERLAY 

Some  minor  additions  to  the  overlay  structures  of  links  8  and  13 
are  required.  Since  the  materials  information  specified  on  MATPZ1  and 
MATPZ2  cards  is  returned  to  the  appropriate  routines  via  named  common 
/MATPZ/,  this  nrme  must  be  added  to  the  overlay  structure  of  links  8  and 
13  and  may  be  included  in  the  same  segment  as  named  commons  /MATIN/  and 
/MATOUT/,  the  standard  materials  named  common.  Also,  in  link  13,  new 
named  common  /ISAVE/  should  be  added  to  the  same  overlay  segment  as 
/SDR2X7/.  This  new  common  allows  the  transfer  of  information  between 
subroutines  SDR2E  and  STPAX3  or  STRAX3. 

Finally,  the  increased  lengths  of  subroutines  in  link  13  (STPAX1 , 
STPAX2,  STPAX3,  STRAX1,  STRAX2,  and  STRAX3)  necessitate  checking  open 
core  named  commons  /SDRA2/  and  /SDR2XX/  for  proper  positioning.  This 
consideration  is  automatically  accounted  for  in  link  8  by  EMG. 
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NOTES 


A  few  notes  concerning  piezoelectric  problems  with  NASTRAN  will 
probably  be  useful. 

1.  In  order  to  use  piezoelectric  materials,  SYSTEM(78)  must  be  set 
to  1  or  2  on  the  NASTRAN  card.  (The  default  value  is  0.)  A  value  of  one 
indicates  electrical-structural  coupling,  a  value  of  2  allows  the  use  of 
piezoelectric  materials  but  does  not  take  into  account  any  electrical 
effects.  The  latter  case  requires  that  the  degrees-of-f reedom 
corresponding  to  the  electric  potential  be  constrained. 

2.  The  electric  potential  at  each  ring  is  considered  to  be  degree- 
of-freedom  A.  Degrees-of-f reedom  5  and  6  always  have  zero  stiffness  and 
must  be  removed  from  the  problem  with  SPCAX  or  RINGAX  cards.  Electroded 
surfaces  (surfaces  of  constant  potential)  may  be  specified  with  MPCAX 
cards . 

3.  Only  TRAPAX  and  TRIAAX  elements  may  reference  piezoelectric 
material  cards  MATPZ1  and  MATPZ2 . 

A.  Standard  material  cards  MAT1  and  MAT3  are  allowed  in  problems 
which  also  contain  piezoelectric  materials. 

5.  The  SE  and  d  values  on  MATPZ1  cards  will  be  multiplied  by  10~12 
by  NASTRAN.  Also  the  value  of  is  fixed  in  NASTRAN  as  8.85A  x  10"-^ 
farad/meter . 

6.  As  may  be  seen  in  Equation  (3),  the  lower  right-hand  portion  of 
the  stiffness  matrix  is  negative-definite.  This  situation  does  not 
affect  NASTRAN  execution  except  that  grid  point  singularity  warning 
messages  are  issued  for  all  unconstrained  electric  potentials. 

7.  To  specify  surface  charge  loads,  the  first  ring  identification 
number  on  the  PRESAX  card  (field  A)  must  be  negative.  This  format 
change  will  allow  NASTRAN  to  distinguish  between  electrical  charges  and 
structural  pressures  within  a  piezoelectric  run.  This  change  is  allowed 
only  when  SYSTEM(78)  =  1. 

8.  Lumped  mass  and  consistent  mass  are  available  for  TRAPAX  and 
TRIAAX  elements.  The  mass  associated  with  the  electric  potential  degree- 
of-freedom  is  zero.  Therefore,  if  a  normal  modes  analysis  by  GIVENS 
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method  is  run,  all  unconstrained  electric  potentials  must  appear  on  OMIT 
cards . 

9.  If  a  structural  damping  coefficient  is  specified  on  a  MATl’Z 1  or 
MATPZ2  card  in  a  dynamics  problem,  the  terms  of  the  resulting  structural 
damping  matrix  corresponding  to  electric  potentials  will  be  zero.  The 
uniform  structural  damping  parameter  G  in  direct  frequency  response 
problems  should  not  be  used,  since  its  use  will  result  in  structural 
damping  terms  corresponding  to  the  electric  potentials. 

10.  Prior  to  the  present  work,  NAS TRAN  could  not  handle  stresses  or 
forces,  whether  real  or  complex,  in  axisymmetric  (AX1C)  dynamics  problems 
NAS TRAN  can  now  handle  all  such  uses  for  the  TRAPAX  and  TRIAAX  finite 
elements . 

11.  Material  properties  specified  on  MATPZ1  cards  are  transformed  by 
NAS TRAN  from  the  standard  1,  2,  3  material  directions  to  the  R,  Z,  0 
directions.  Also,  the  transformation  required  due  to  a  switch  in  the 
order  of  the  R-Z  and  Z-Q  shears  between  conventional  specifications  and 
NASTRAN  is  performed  for  MATPZ1  properties.  Material  properties  on 
MATPZ2  cards  are  used  by  NASTRAN  as  they  appear  on  the  card.  Therefore, 
any  required  transformation  must  be  performed  by  the  user. 


SAMPLE  PROBLEMS 

Two  sample  problems  were  ran  to  check  results  from  TRIAAX  and  TRAPAX 
piezoelectric  finite  elements  against  experimental  and  other  calculated 
results . 

Problem  1 

The  natural  frequencies  of  a  radially  polarized  PZT-4  hollow  cylindo 
were  computed.  Figure  5  shows  the  finite  element  mesli  and  boundary 
conditions  of  the  problem  described  by  Allik.^  Results  from  the  present 
work  were  compared  with  the  results  from  the  mesh  in  Figure  5. 

For  the  present  analysis  of  the  hollow  cylinder,  four  finite  element 
meshes  were  used.  Mesh  1,  shown  in  Figure  6,  consists  of  3b  TRIAAX 
elements  with  45  nodes.  One  other  mesh  with  TRIAAX  elements  was 
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analyzed,  quadrupling  the  number  of  elements  of  the  first  mesh.  Mesh  2 
consists  of  224  TRIAAX  elements  and  145  nodes.  One  analysis  was  done 
using  TRAPAX  elements.  Mesh  3  consists  of  112  TRAPAX  elements  and  145 
nodes . 

For  each  mesh  two  analyses  were  done,  the  first  considering  the 
effect  of  piezoelectric  coupling  and  the  second  neglecting  piezoelectric 
coupling.  The  results  of  all  the  analyses  are  presented  in  Tables  1  and 
2.  Allik's  results  were  estimated  from  a  graph  in  that  article.^  Those 
results  were  obtained  using  quadratic  isoparametric  triangular  finite 
elements  and  the  MARTSAM  program. ^ 


TABLE  1  -  NATURAL  FREQU!  CIES  OF  HOI, LOW  CYLINDER 
INCLUDING  PIEZOELECTRIC  COUPLING 


Natural  Frequencies  (Rad/Sec  x  10 

’) 

Mode 

Mesh  from  Ref.  (1) 

Mesh  1 

Mesh  2 

Mesh  3 

1  (R) 

2.88 

2.91 

2.90 

2.90 

2  (R) 

3.25 

3.  37 

3.21 

3.23 

l  (D 

5.13 

5.14 

5.12 

5.08 

3  (R) 

7.25 

7.90 

7.05 

7.20 

4  (R) 

14.0 

1  3.9 

13.1 

13.5 

r1 

15.0 

14.9 

. 

14.3 

14.4 

(R)  denotes  radial  mode 

(L)  denotes  longitudinal  mode 


Problem  2 

The  natural  frequencies  of  an  axial Iv  polarized  1V.T-4  piezoelectric 
disk  were  computed.  This  disk  is  shown  in  Figure  7.  Both  the  upper 
and  lower  electroded  surfaces  are  grounded.  The  results  ol  the  present 
work  were  compared  to  calculated  and  experimental  results  of  the  same 
disk  provided  by  the  Naval  Underwater  Systems  Center  (NUSC)  . 

The  mesh  for  the  present  analysis  of  the  piezoelectric  disk  contains 
205  nodes  and  320  TRIAAX  elements.  The  NUSC  model  contained  83  nodes  and 
28  quadratic  isoparametric  triangular  elements  and  was  run  with  the 
MARTSAM  program.  The  results  of  the  analyses  are  presented  in  Table  3. 


TABLE  2  -  NATURAL  FREQUENCIES  OF  HOLLOW  CYLINDER 
NEGLECTING  PIEZOELECTRIC  COUPLING 


Mode 

Natural  Frequencies  (Rad/Sec  x  10‘ 

-') 

Mesh  from  Ref.  (1) 

Mesh  1 

Mesh  2 

Mesh  3 

1  (R) 

2.88 

2.90 

2.90 

2.90 

2  (R) 

3.25 

3.37 

3.21 

3.16 

l  (L) 

5.13 

5.09 

5.08 

5.08 

3  (R) 

6.50 

7.58 

6.69 

6.34 

4  (R) 

11.9 

13.1 

11.9 

11.4 

(R)  denotes  radial  mode 

(L)  denotes  longitudinal  mode 


TABLE  3  -  NATURAL  FREQUENCIES  OF  PIEZOELECTRIC  DISK 


Mode 

Natural  Frequencies  (cp 

a) 

NUSC  Mesh 
Exper imental 

NUSC  Mesh 
Calculated 

NASTRAN 

Mesh 

1 

22042 

23298 

24323 

2 

59805 

61114 

3 

103048 

104689 
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